Acquisition: San Francisco (USA), ALOS PALSAR, L-band
Path to images: /projects/data/polsar/ALOS-P1_1A-ORBITALPSRP202350750/
SLC (single-look complex) real & imag parts files:
Original Image size (px in azimuth x range): 18432 x 1248
Tips:
%matplotlib widget
# import useful libraries, functions, and modules
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import uniform_filter
Step 1 : Load data
# path to data
path = '/projects/data/polsar/ALOS-P1_1__A-ORBIT__ALPSRP202350750/'
img_shape = (18432, 1248)
def read_slc(pol_name):
slc = np.fromfile(path + 'i_%s.bin' % pol_name, dtype=np.float32) + 1j*np.fromfile(path + 'q_%s.bin' % pol_name, dtype=np.float32)
return slc.reshape(img_shape) / 1e6
# Read images
slchh = read_slc('HH')
slchv = read_slc('HV')
slcvh = read_slc('VH')
slcvv = read_slc('VV')
# Print image size
slchh.shape
(18432, 1248)
# Crop images in azimuth
crop_az = (7000, 12000)
slchh = slchh[crop_az[0]:crop_az[1], :]
slchv = slchv[crop_az[0]:crop_az[1], :]
slcvh = slcvh[crop_az[0]:crop_az[1], :]
slcvv = slcvv[crop_az[0]:crop_az[1], :]
# Print image size
slchh.shape
(5000, 1248)
Visualize the HH amplitude image before & after multilook and compare intensity PDF distributions
Amplitude corresponds to the magnitude of the complex SLC (np.abs(slc)):
Intensity corresponds to the squarred amplitude (np.abs(slc)**2):
Tips:
looksa = 9
looksr = 5
# Multi looking in HH
int_hh = np.abs(slchh)**2
amp_hh_ml = np.sqrt(uniform_filter(int_hh, [looksa, looksr]))
del int_hh
# plot
plt.figure(figsize = (10,10))
ax = plt.subplot(2,1,1)
plt.imshow(np.transpose(np.abs(slchh)), vmin = 0, vmax = 3*np.mean(np.abs(slchh)), cmap = 'gray', aspect = 'auto')
ax.set_title('Single look HH')
plt.tight_layout()
ax = plt.subplot(2,1,2)
plt.imshow(np.transpose(amp_hh_ml), vmin = 0, vmax = 3*np.mean(amp_hh_ml), cmap = 'gray', aspect = 'auto')
ax.set_title('Multi look HH')
plt.tight_layout()
Plot the PDF (histogram) of the intensity at HH before and after the Multilook
Tips:
# plot
plt.figure(figsize = (8,8))
ax = plt.subplot(2,1,1)
ax.hist((np.abs(slchh[:,:200])**2).ravel(), 100)
ax.set_title("PDF original SLC")
ax = plt.subplot(2,1,2)
ax.hist((amp_hh_ml[:,:200]**2).ravel(), 100)
ax.set_title("PDF After Multilook")
plt.tight_layout()
Create, visualize and compare (where are the same? where are different?) RGB composite images using
Tips:
#### Lexicographic and Pauli representations
# Lexicographic
amphh = np.sqrt(uniform_filter(np.abs(slchh)**2, [looksa, looksr]))
ampvv = np.sqrt(uniform_filter(np.abs(slcvv)**2, [looksa, looksr]))
amphv = np.sqrt(uniform_filter(np.abs(slchv)**2, [looksa, looksr]))
# Initialize RGB array (Npx_rg, Npx_az, 3)
dimaz = amphh.shape[0]
dimrg = amphh.shape[1]
rgb_lex = np.zeros((dimrg, dimaz, 3), 'float32')
rgb_lex[:,:,0] = np.clip(np.transpose(amphh), 0, 2.5*np.mean(amphh))
rgb_lex[:,:,1] = np.clip(np.transpose(amphv), 0, 2.5*np.mean(amphv))
rgb_lex[:,:,2] = np.clip(np.transpose(ampvv), 0, 2.5*np.mean(ampvv))
rgb_lex[:,:,0] = rgb_lex[:,:,0] / np.max(rgb_lex[:,:,0])
rgb_lex[:,:,1] = rgb_lex[:,:,1] / np.max(rgb_lex[:,:,1])
rgb_lex[:,:,2] = rgb_lex[:,:,2] / np.max(rgb_lex[:,:,2])
# delete variables
del amphh, amphv, ampvv
# Pauli amplitudes
pauli1 = np.abs(slchh+slcvv)**2
pauli1 = np.sqrt(uniform_filter(pauli1, [looksa, looksr]))
pauli2 = np.abs(slchh-slcvv)**2
pauli2 = np.sqrt(uniform_filter(pauli2, [looksa, looksr]))
pauli3 = np.abs(2*slchv)**2
pauli3 = np.sqrt(uniform_filter(pauli3, [looksa, looksr]))
# initialize array
rgb_pauli = np.zeros((dimrg, dimaz, 3), 'float32')
rgb_pauli[:,:,0] = np.clip(np.transpose(pauli2), 0, 2.5*np.mean(pauli2))
rgb_pauli[:,:,1] = np.clip(np.transpose(pauli3), 0, 2.5*np.mean(pauli3))
rgb_pauli[:,:,2] = np.clip(np.transpose(pauli1), 0, 2.5*np.mean(pauli1))
rgb_pauli[:,:,0] = rgb_pauli[:,:,0] / np.max(rgb_pauli[:,:,0])
rgb_pauli[:,:,1] = rgb_pauli[:,:,1] / np.max(rgb_pauli[:,:,1])
rgb_pauli[:,:,2] = rgb_pauli[:,:,2] / np.max(rgb_pauli[:,:,2])
# plots
plt.figure(figsize = (12,12))
## lexicographic
ax = plt.subplot(2,1,1)
plt.imshow(rgb_lex, aspect = 'auto')
ax.set_title('RGB lexicographic basis')
plt.tight_layout()
## Pauli
ax = plt.subplot(2,1,2)
plt.imshow(rgb_pauli, aspect = 'auto')
ax.set_title('RGB Pauli basis')
plt.tight_layout()
Computation of the alpha1 angle.
del slchh, slcvv, slchv
# Compute the alpha angle
# --- calculate the length of the Pauli vector
pauli_l = np.sqrt(np.abs(pauli1)**2 + np.abs(pauli2)**2 + np.abs(pauli3)**2)
# --- compute alpha angle [rad]
alpha = np.arccos(np.abs(pauli1) / pauli_l)
# plots
plt.figure(figsize=(12,12))
plt.subplot(2,1,1)
plt.imshow(rgb_pauli, aspect='auto')
plt.colorbar() # dummy colorbar to align images
plt.subplot(2,1,2)
plt.imshow(np.transpose(alpha) * 180/np.pi , cmap = 'jet', vmin = 0, vmax = 90, aspect = 'auto', interpolation = 'nearest')
plt.colorbar()
plt.tight_layout()